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Chapter  1 


Introduction 

1.1  Position  and  Velocity  Estimation 

Wireless  location  finding  has  emerged  as  an  essential  public  safety  feature  in  cellular  networks 
due  to  a  mandate  issued  by  the  FCC  in  1996  [1].  The  E91 1  mandate  requires  wireless  service 
providers  to  deliver  the  accurate  location  of  a  91 1  caller  to  emergency  responders.  In  addition  to 
its  immediate  EMS  and  military  applications,  there  exists  much  potential  in  the  commercial 
sector  as  well  [1][2].  As  the  use  of  location  sensitive  applications  becomes  more  prevalent,  the 
demand  for  reliable,  cost-effective  passive  location  systems  will  subsequently  grow  [1][2]. 

A  commonly  used  method  for  estimating  the  position  of  an  emitter  within  a  wireless  sensor 
network  requires  the  use  of  time-difference-of-arrival  (TDOA)  measurements  [3] [4],  The  locus  of 
points  where  the  difference  in  distance  to  two  sensors  is  proportional  to  the  TDOA  estimate 
obtained  from  the  sensors  is  a  hyperbola  on  which  the  emitter  lies  (Eig.  1.1).  If  the  elevation  of 
the  emitter  is  known,  at  least  3  sensors  are  required  to  determine  its  position  [5].  Including 
additional  information  from  redundant  sensors  can  help  to  improve  the  accuracy  of  the 
estimation.  Similarly,  the  velocity  of  an  emitter  can  be  determined  using  frequency-dijf erence-of- 
arrival  (FDOA)  information. 


Eigure  1.1  Passive  location  of  emitter  using  TDOA  estimates 


1.2  The  Data  Compression  Issue 

In  order  to  estimate  TDOA  or  EDOA  parameters,  it  is  necessary  for  information  to  be  shared 
between  sensors.  Sharing  data  among  sensors  can  introduce  significant  network  delays  that  may 
affect  the  performance  of  the  location  application.  To  minimize  network  latency  and  reduce 
power  consumption  by,  possibly,  mobile  sensors,  the  data  collected  by  each  sensor  should  be 
compressed  before  it  is  transmitted  to  a  neighboring  sensor  (Eig  1.2a)  or  a  central  fusion  center 
(Eig  1.2b)  for  processing. 
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Figure  1.2  (a)  distributed  estimation  (b)  eentralized  estimation 

In  addition  to  noise  introduced  at  each  receiver,  the  combined  decimation  and  quantization  of  the 
signal  information  will  also  contribute  noise  to  TDOA  estimate.  Unlike  traditional  perceptual 
encoding  schemes,  the  fidelity  criterion  of  this  particular  application  is  the  impact  on  the 
TDOA/FDOA  parameter  estimates. 


1.3  Scope  of  Work 

This  thesis  will  outline  the  issues  related  with  compressing  data  for  the  purpose  of  geo-location 
and  will  develop  a  tool  which  will  aid  in  the  design  of  compression  algorithms  for  digital  linearly 
modulated  signals  of  the  form; 


v(t)=  Y.hgit-nT)  (1.1) 

«=-oo 

Where  {/,,}  is  a  sequence  of  symbols,  gj{t)  is  a  transmission  filter,  and  v{t)  is  the  continuous 
time  complex  envelope  of  the  modulated  signal.  Examples  of  digital  modulation  techniques 
which  are  of  the  form  in  equation  1 . 1  include  pulse  amplitude  modulation  (PAM),  phase  shift 
keying  (PSK),  and  quadrature  amplitude  modulation  (QAM). 

The  developments  in  the  remainder  of  this  thesis  do  not  include  the  effects  of  multipath  or 
fading.  Additionally,  it  is  assumed  that  the  sensors  and  fusion  center  have  clocks  which  are 
perfectly  synchronized,  so  the  effects  of  clock  jitter  are  also  not  included. 
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Chapter  2 


Literature  Review 

2.1  Data  Compression  Overview 

Data  Compression  or  Source  Coding  is  the  process  of  removing  unwanted  redundancy  from  a 
signal  while  satisfying  a  fidelity  criterion  to  reduce  transmission  bandwidth  or  space  required  for 
storage  [6].  Some  examples  of  popular  perceptual  [6]  encoding  schemes  include  the  MPEG 
[6][7]  and  JPEG  [6][7]  standards. 

A  common  lossy  data  compression  algorithm  structure  consists  of  several  stages  as  depicted  in 
Eig.  2.1.  Eirst,  the  input  is  transformed  to  a  domain  suitable  for  encoding.  The  samples  of  the 
transformed  signal  are  then  selectively  mapped  to  a  smaller  set  of  symbols,  which  require  fewer 
bits  to  represent,  through  the  process  of  quantization.  The  compressed  data  is  then  transmitted 
over  a  channel  to  a  decoder  which  decodes  the  quantized  sequence  and  restores  the  information 
to  its  original  domain. 


Eigure  2.1  A  common  lossy  data  compression  algorithm  structure 

The  majority  of  existing  lossy  compression  schemes  are  designed  for  multimedia  applications 
and,  consequentially,  use  perceptual  based  measures  to  determine  loss  of  fidelity.  The  focus  of 
this  thesis  will  be  on  the  effects  of  lossy  compression  on  TDOA/EDOA  estimation. 


2.2  Transform  Coding 

It  is  common  in  data  compression  schemes  to  transform  a  signal  into  a  domain  that  is  more 
suitable  for  analysis  or  encoding.  Eor  example,  in  subband  coding,  a  signal  is  decomposed  into  a 
number  of  frequency  subbands  (Eig.  2.2),  each  of  which  can  be  encoded  independently  (Eig  2.2). 
There  exist  many  techniques  for  performing  subband  decomposition  including  quadrature -mirror 
filtering  (QME)  and  polyphase  decomposition  [7].  Two  special  cases  of  these  methods  are 
presented;  the  wavelet  packet  transform  [7]  and  the  pseudo-quadrature  mirror  filter  bank  [8]. 
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Figure  2.2  Subband  Decomposition 


The  wavelet  packet  transform  (WPT)  is  a  cascade  of  quadrature -mirror  filter  banks  and  can  take 
on  any  dyadic  tree  structure  (Fig  2.3).  Because  QMF  banks  satisfy  the  perfect  reconstruction 
perfect  reconstruction  conditions,  the  input  to  a  WPT  analysis  filter  bank  can  be  perfectly 
reconstructed  by  a  companion  synthesis  filter  bank.  Aliasing  which  is  introduced  by  overlapping 
from  adjacent  frequency  bands  will  be  completely  cancelled.  Although  the  WPT  allows  for 
perfect  reconstruction,  due  to  its  recursive  nature,  it  is  not  the  most  computationally  efficient 
structure. 


Figure  2.3  Wavelet  Packet  Decomposition 


If  the  perfect  reconstruction  conditions  are  relaxed,  a  more  efficient  structure  can  be 
implemented  such  as  the  pseudo-quadrature  mirror  fdter  bank  (PQF)[8].  The  pseudo-quadrature 
mirror  filter  bank  is  a  cosine-modulated  filter  bank  which  is  designed  to  minimize,  but  not 
eliminate,  aliasing  between  adjacent  frequency  subbands.  Section  2.6  provides  a  more  detailed 
description  of  the  filter  bank  and  a  brief  analysis  of  its  performance.  In  general,  the  order  of  a 
near-perfect  reconstruction  filter  is  much  smaller  than  that  of  a  perfect  reconstruction  filter  with 
the  same  transition  band  slope.  Additionally,  the  number  of  operations  required  by  a  polyphase 
filter  bank  is  less  than  that  of  a  recursive  filter  bank. 


2.3  Quantization 

Quantization  is  the  process  of  mapping  the  outputs  of  an  information  source,  x,  to  a 
representation  sequence,  x,  which  requires  fewer  bits  to  describe.  The  function  which  defines 
the  mapping  is  called  a  quantizer  [9]  [10].  In  scalar  quantization,  each  single-source  output  is 
quantized  into  a  number  of  levels  and  the  levels  are  encoded  into  a  binary  sequence  [9]  [10].  If 
we  are  dealing  with  an  analog  source,  each  output  is  a  real  number.  A  scalar  quantizer,  Q{  )  , 
partitions  91  into  N  disjoint  subsets,  9I„1  <i<N  [11].  Corresponding  to  each  subset,  91,,  is  a 
representation  point,  x,  (Fig.  2.4). 
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Figure  2.4  Scalar  quantization  of  91 

If  the  source  output  at  time  k,  belongs  to  91,,  it  is  then  represented  by  x,  [9]. 

Q{x)  =  X.  for  all  x  e  91,.  (2.1) 


The  quantized  symbol  at  time  k,  Xj^,  is  then  encoded  as  a  binary  sequence  and  transmitted.  Since 
there  are  N  quantization  levels,  the  number  of  bits  required  to  transmit  each  source  output  is 
R  =  \og2N. 


Since  quantization  provides  a  many-to-one  mapping  (i.e.  91,  x,),  the  original  source  symbol, 

Xj.,  can  not  be  obtained  from  its  quantized  version,  The  distance  between  a  source  output  at 
time  k,  x,^,  and  its  quantized  version,  X/^,  can  be  expressed  as  the  function  J(xj.,x,5.).  A 
commonly  used  distortion  measure  is  the  squared-error  [6]  [9]  [1 1][12]; 

d(Xk,x,)  =  ix,-x,y  (2.2) 


If  d{Xj^,Xj^)  is  the  distortion  measure  per  single  source  output,  then  distortion  between  a  sequence 
of  n  outputs,  X” ,  and  the  corresponding  quantized  values,  X"  ,is  the  average  of  the  distortion 
between  the  individual  outputs: 


d{X\X')^\Y.d{x„x,)  (2.3) 

k=\ 


Equation  2.3  assumes  that  the  position  of  the  error  in  the  sequence  is  unimportant  [9].  Since  the 
source  output  is  a  random  process,  t/(X",X'')  is  a  random  variable.  The  distortion  for  the  source 
is  defined  as  the  expected  value  of  the  random  variable  (i(X",X"). 


D  =  EidiX'X)-]  =  \Y.E[d{x,X)'\  =  Eidix.X)^  (2.4) 

k=\ 

Where  the  final  equality  in  equation  2.4  holds  under  the  assumption  that  the  source  output  is  a 
stationary  random  process  [12]. 

The  distortion,  D,  can  be  expressed  as  a  function  of  the  rate,  [9][1 1].  In  general,  D  is 
monotonically  decreasing  in  i?  [9][10][1 1][12].  A  well  known  result  in  rate-distortion  theory 
[10]  [12]  is  that  the  distortion-rate  function,  £)(i?),  of  a  discrete-time,  continuous-amplitude, 
memoryless  source  with  zero-mean  and  finite  variance,  aj',  (with  respect  to  the  mean-squared- 
error  distortion  measure)  is  upper  bounded  by: 
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(2.5) 


DiR)  <  0-/2-''' 

Where  equality  holds  when  the  source  outputs  are  Gaussian  with  variance  o-^'  [12],  The  limit  in 
2.5  can  be  approached  asymptotically  using  very  complex  encoders  and  decoders  [9]. 


The  simplest  scalar  quantizer  is  the  uniform  quantizer,  which  has  regions,  91 .,  of  equal  length.  A, 
(with  the  exception  of  9ij  and  91^).  Uniform  quantization  is  optimal  when  the  probability 
density  function  (PDF)  of  the  source  output  is  uniform.  If  the  PDF  of  the  source  output  is  known 
and  is  not  uniform,  an  optimal  quantizer  can  be  designed  through  use  of  an  algorithm  which 
iteratively  determines  the  boundaries  of  the  partition  regions  [13]. 

For  the  scenario  of  a  wireless  sensor  network,  we  use  uniform  scalar  quantization  due  to  its  ease 
of  implementation  and  our  lack  of  knowledge  of  source  statistics.  For  a  finite  amplitude  signal 
in  additive  noise,  the  partition  length.  A,  is  given  by  [19]; 


A  = 


r^b-i  _  j_ 
^  2 


(2.6) 


Where  A  is  the  peak  amplitude  of  the  signal,  cr'  is  the  variance  of  additive  noise,  and  b  is  the 
number  of  bits  allocated  for  quantization. 


2.4  Bit  Allocation 

In  a  transform  coding  scheme,  it  is  often  required  to  encode  segments  of  a  signal  independently, 
for  example,  in  subband  encoding,  each  frequency  subband  is  separately  encoded.  Given  a 
single  rate  constraint  and  a  set  of  subbands  to  be  individually  quantized,  the  problem  of 
determining  an  appropriate  way  to  allocate  bits  arises. 

As  mentioned  in  the  previous  section,  the  distortion  introduced  by  quantization  can  be  expressed 
as  a  function  of  the  rate,  D{R).  For  an  M-channel  subband  encoder,  each  subband  will  have  its 
own  rate-distortion  function,  D^{R),\<  i<  M .  When  the  distortion  measure  is  additive,  the  total 
distortion  can  be  computed  as  the  sum  of  the  distortion  of  the  individual  subbands.  Given  a  set 
of  independent  rate-distortion  functions,  D.{R),  the  task  of  finding  a  bit  allocation  vector, 

B  =  \b.  >  0|/  e  1,2,-  •  which  minimizes  the  total  distortion,  using  less  than  R  bits,  is  defined 
as: 


,<i?j  (2.6) 

When  a  signal  is  stationary  and  its  statistics  are  known  a  priori,  a  model  can  be  used  to  determine 
a  bit  allocation  vector,  B,  which  will  be  optimal  on  average.  If  the  signal  is  highly  non-stationary 


M 


subjec 


t  to; 


M 


/=! 
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or  its  statistics  are  not  known,  an  operational  algorithm  ean  be  used  to  ealeulate  an  optimal  bit 
alloeation  for  each  bloek  whieh  is  quantized.  Determining  the  bit  alloeation  requires  solving  the 
eonstrained  optimization  problem  in  equation  2.6,  whieh  may  not  be  eomputationally  feasible  for 
a  wireless  sensor  to  perform  on  eaeh  quantization  bloek.  An  effieient  Lagrange  multiplier  based 
teehnique  [14]  has  been  implemented  to  solve  the  integer  optimization  problem  off-line  (see 
Appendix  A).  It  is  assumed  that  a  eentral  fusion  eenter  will  observe  a  portion  on  an  intereepted 
signal  and  ealeulate  a  bit  alloeation  veetor,  B,  whieh  will  then  be  sent  to  the  sensors.  Performing 
the  bit  alloeation  at  the  fusion  center  will  reduee  the  resourees  required  by  eaeh  sensor  without 
signifieantly  affeeting  the  overall  geo-loeation  task. 

If  the  rate-distortion  funetions  are  not  eonvex,  or  do  not  eontain  a  point  whieh  meets  the  required 
budget,  then  the  Lagrange  optimization  method  will  not  provide  an  optimal  solution  [6].  When 
the  Lagrange  multiplier  method  is  inadequate,  the  alloeation  ean  be  formulated  as  a  deterministie 
dynamie  programming  problem.  Although  dynamie  programming  ean  provide  a  better  solution, 
it  requires  signifieantly  more  resourees  to  implement  and  may  not  be  feasible  in  real-time 
applieations. 


2.5  TDOA  and  FDOA  Parameter  Estimation 

In  order  to  develop  a  distortion  measure  whieh  is  appropriate  for  geo-loeation,  it  is  important  to 
examine  the  TDOA  and  FDOA  estimation  proeess.  The  TDOA  parameter  is  determined  by 
observing  data  reeeived  by  a  pair  of  physieally  separated  sensors  (Fig  2.5a)  and  searehing  for  the 
time  offset  of  a  hypothesized  eommon  signal.  Likewise,  the  FDOA  parameter  ean  be  found  by 
determining  the  frequeney  offset  between  the  signals  reeeived  at  the  sensors  (Fig.  2.5b).  A 
eommon  way  to  perform  the  estimation  proeess  is  through  the  use  of  eross-eorrelation  [15]. 


sensor  1 


o 


time  difference  of  arrival 


sensor  1 


o 


sensor  2 


emitter 


* 


frequency  difference  of  arrival 


Q. 


sensor  2 


emitter 


(a)  (b) 

Figure  2.5  (a)  stationary  emitter  (b)  mobile  emitter  in  wireless  sensor  network 


Let  Xj(t) denote  the  signal  observed  by  sensor  1  and  X2(0the  signal  observed  by  sensor  2.  Let 
both  signals  eontain  a  eommon  eomponent  5(t)and  a  white  Gaussian  noise  eomponent,  n.{t) , 
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which  is  assumed  to  be  independent  at  eaeh  sensor. 

Xj(t)  =  +  n^{t)  (2.7) 

x^it)  =  sit-D)e^^’^^‘  +  n^(t)  (2.8) 

The  complex  ambiguity  function  (CAP)  [15]  [4]  is  defined  as: 

A{T,f)  =  \  \ft)xft+  T)e~^^’^‘dt  (2.9) 

Notice  that  the  CAP  is  the  eross-correlation  of  the  signals  observed  by  the  two  sensors, 
modulated  by  a  term  f.  The  time  lag,  t  ,  and  frequency  shift,/,  parameters  are  to  be  searched 
simultaneously  for  values  that  maximize  \A(t,  /)| .  The  combined  TDOA/PDOA  estimation 
process  can  be  expressed  as: 


max|y4(r,/)|  (2.10) 


Assuming  that  there  are  no  hidden  periodicities  in  the  signals  Xj(t)  and  X2it),  and  the  parameters 
r  and / remain  constant  within  the  integration  time  f;  the  ambiguity  function  will  peak  as 
T  ^  D  and  /  f\-  f 2^  creating  a  unique  “correlation  lobe”  (Pig.  2.6). 


Complex  Ambiquity  Function 


Pigure  2.6  Complex  Ambiguity  Punction 
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Given  that  s(t)  is  aperiodic,  the  CAP  should  possess  a  unique,  clearly  defined  maximum.  When 
the  signals  x,(0  contain  additive  noise  n^t),  the  peak  location  of  the  CAP  will  be  perturbed.  It 
has  been  shown  [4][15][16]  that  the  TDOA  and  PDOA  estimates  are  unbiased  and  have  a 
variance  that  achieves  the  Cramer-Rao  bound  when  any  of  several  reasonable  techniques  is  used 
to  locate  the  apparent  peak.  When  the  additive  noise,  n.(t),  is  white  and  Gaussian,  the  Cramer- 
Rao  lower  bounds  for  TDOA  and  PDOA  estimation  are  given  as  [15][16]; 


^TDOA  — 


^FDOA  — 


1 

1 

where  (3  =  271 

p^BTy 

1 

1 

where  T  =27i 

T.  fBTy 

e 

J  f\S{f)fdf 

.  J  |5(/)r<r . 

s(t)^dt 


(2.11) 


(2.12) 


I  |2 

Where  |5'(/)|  is  the  power  spectrum  density  of  the  source  signal  s(t),B  is  the  noise  bandwidth 
at  the  receiver,  T  is  the  integration  time,  fd  is  termed  rms  radian  frequency,  is  termed  rms 
integration  time,  and  y  is  the  effective  input  signal-to-noise  ratio  (SNR)[4][15]; 


r 


2 


1  1 

—  +  — 

.n  Yi 


+ 


1 

YxYi- 


(2.13) 


Where  y^  is  the  SNR  at  sensor  1  and  is  the  SNR  at  sensor  2.  In  order  to  accurately  perform 
the  estimation,  the  integration  time-bandwidth  product  must  be  large,  BT  » 1 .  A  small  effective 
SNR,  y ,  will  degrade  the  accuracy  of  the  estimate;  however,  the  accuracy  can  be  further 
improved  by  increasing  the  receiver  bandwidth,  B,  or  observation  time,  T. 

In  order  to  improve  the  asymptotic  performance  of  an  unbiased  estimator,  we  wish  to  minimize 
the  Cramer-Rao  lower  bound  or,  equivalently,  maximize  its  reciprocal,  the  Pisher  Information,  J 
[17]. 


(7^ TDOA  >  /  *  (2.14) 

Lossy  compression  of  the  observed  signal,  x.(t),  will  change  its  probability  density  function 
(PDP),  which  will,  consequently,  affect  the  value  of  the  Pisher  Information  (PI).  Let  J(x;  r) 
denote  the  PI  of  the  observed  signal  with  respect  to  the  TDOA  parameter,  r,  and  J(x;  r)  denote 
the  PI  of  the  compressed  signal.  Our  goal  is  to  compress  the  observed  signal  in  such  a  way  that 
the  difference  in  Pisher  Information  is  minimized; 

min{J(x;r)- J(x;r)}  (2.14) 


Since  J(x;  r)  >  J(x;  r),  for  any  compression  scheme,  the  goal  is  [19]: 
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max{  J(x;  r)} 


(2.15) 


a: 

The  Fisher  Information  for  TDOA  estimation  when  the  two  observed  signals  contain  additive 
white  Gausian  noise  can  be  expressed  as  [15][16]; 


fsifYdf 


J{x-,t)  =  BTy- 


J  |5(/)f 


df 


(2.16) 


When  the  SNR  at  the  two  sensors  is  large,  i.e.  y  »  h/j  »  1,  the  third  term  in  equation  2.13  is 
negligible  and  the  effective  SNR,  /,  is  approximately: 


n  +  Vi 


(2.17) 


The  SNR  for  the  sensor  is  defined  as: 


J  Wtdf 

r,=- — ^ —  (2.18) 

Where  a]  is  the  variance  of  the  zero-mean  additive  white  Gaussian  noise  at  the  sensor.  Using 
equation  2.18,  the  effective  SNR,  /,  can  be  rewritten  as: 

2 1  \Sif)fdf 

r=^2 — ^  (2-19) 

Now,  substituting  the  expression  for  the  effective  SNR  (Eq.  2.19)  into  equation  2.16,  the  Fisher 
Information  can  then  be  expressed  as: 


J  f-\S(f)fdf 

J{x;  t)  =  2BT^ — ^ ^ -  (2.20) 

(Tj  +  a  2 

Focusing  on  the  compression  task  at  sensor  i,  we  wish  to  compress  the  intercepted  signal,  x^{f), 
to  maximize  equation  2.20.  For  a  fixed  time-bandwidth  product,  BT,  the  term  which  we  wish  to 
maximize  at  sensor  i  is  proportional  to: 


y]f\S{f)fdf  (2.21) 
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This  approach  has  been  used  in  previous  works  [18][19][20][22]  in  whieh  the  souree  signal,  ^(t), 
is  a  frequeney  modulated  (FM)  signal.  Their  works  indieate  that  the  FI-based  approaeh  provides 
better  results  than  using  standard  MSB  distortion  eriteria  when  eompressing  signals  for  TDOA 
and  FDOA  estimation.  For  the  remainder  of  this  thesis,  we  will  follow  the  work  of  [19]  [21]  [22] 
in  developing  a  model  for  the  Fisher  Information  whieh  eharaeterizes  the  effeets  of  lossy 
eompression  on  the  TDOA  estimation  task. 

As  diseussed  in  seetion  2.3,  quantization  introduees  noise.  The  signal,  x,[A:],  resulting  from 
quantization  of  samples  of  the  intereepted  signal,  x.{t),  ean  be  expressed  as; 

x{k]  =  s[k]  +  n,[k]  +  e[k]  (2.22) 

Where  s\k^  is  the  sampled  souree  signal  at  time  k,  n,[A:]  is  the  sensor  noise,  and  e[A:]  is  the 
quantization  noise.  Following  the  work  of  [19][22],  we  will  model  the  quantization  noise,  e,  as 
white,  uniformly  distributed,  zero  mean,  and  independent  of  the  sensor  noise,  n .,  when  the 
samples  are  quantized  using  multiple  bits  [19]  [22].  For  multi-bit  quantization,  if 
var{e}  «  var{n,  },  then  the  sum  of  the  random  variables,  e  and  n,  ,  will  result  in  a  random 
variable  that  is  approximately  Gaussian,  with  a  varianee  that  is  the  sum  of  the  varianees  of  the 
two  random  variables  [19].  Ineluding  the  effeets  of  quantization  noise  (for  the  multi-bit  ease) 
into  equation  2.21,  we  ean  then  express  the  quantity  whieh  we  are  trying  to  maximize  as: 


(2.23) 


Where  is  the  varianee  of  the  quantization  noise,  whieh  we  model  by  [19]; 

(2.24) 

Where  cr^  is  the  varianee  of  the  intereepted  signal,  and  b  is  the  number  of  bits  used  for 
quantization.  As  done  in  [19]  [21]  [22],  we  will  use  the  Diserete  Fourier  Transform  (DFT)  of  the 
intereepted  signal,  x^{t),  to  evaluate  the  Fisher  Information  for  the  multi-bit  quantization  ease 
[19]  (the  superseript  “m”  in  signifies  multi-bit). 


2;r^ 

^ -  (2.25) 

XcJ,  +q„ 

Where  A.[n],-%  <«<%-!  are  the  DFT  eoeffieients  of  the  samples  of  the  signal  X;(t) ;  erf  is  the 
varianee  of  the  sensor  noise,  n,  (0  (scaled  by  Aas  a  result  of  the  DFT  [19]);  and  ql  is  the 
variance  of  the  quantization  noise  for  the  sample.  Notiee  that  the  Fisher  Information  (Eq. 
2.20)  is  dependent  on  the  power  spectrum  density  of  the  source  signal,  s{t) .  For  the  FI 
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evaluation  (Eq.  2.25)  we  are  using  the  DFT  of  the  observed  signal,  which  contains  both 
the  source  signal,  s{t),  and  the  sensor  noise,  therefore,  the  expression  which  we  are 

evaluating  (eq.  2.25)  is  approximately  the  Fisher  Information. 

Since  the  Fisher  Information  is  proportional  to  the  integral  of  the  quadratically  weighted  PSD 
(Eq,  2.21),  the  choice  of  the  DFT  as  the  transform  for  analysis  allows  us  examine  the 
contribution  of  frequency  components  to  the  overall  FI;  therefore,  we  can  selectively  quantize 
frequency  subbands  to  maximize  the  total  FI. 


For  the  case  of  quantization  to  a  single  bit,  the  assumption  that  the  variance  of  the  sensor  noise  is 
large  than  the  variance  of  the  quantization  noise  no  longer  holds.  We  use  the  numerically 
computed  FI  given  in  [19]  for  quantization  using  a  single  bit: 
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Im^{X.[n]}exp 


Re^{X.[n]}exp 
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Re^RW} 
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W{X.[n]} 


'  l-erf\RQ{X,[n]}/i^a)) 

_ 1 _ 

1  -  erf  (Im{X.[n] }  /(ViVo-))_ 
(2.26) 


Where  [«],-%  <«<%-!  are  the  DFT  coefficients  of  the  samples  of  the  signal  X;(t),  and  the 
superscript  “1”  in  signifies  single-bit  quantization. 

In  [21]  the  authors  use  the  Discrete  Fourier  Transform  to  evaluate  the  Fisher  Information  of  an 
intercepted  signal,  while  using  the  Wavelet  Packet  Transform  to  transform  the  data  for  encoding. 
Similar  to  their  approach,  we  will  use  the  DFT  to  evaluate  the  FI,  and  employ  the  Pseudo-QMF 
bank  to  transform  the  data  for  encoding. 


2.6  Pseudo-QMF  Filter  Bank 

The  pseudo-quadrature  mirror  filter  (PQF)  bank  [8]  is  a  cosine-modulated  filter  bank  which  is 
designed  to  cancel  aliasing  between  adjacent  frequency  subbands.  Aliasing  introduced  by  non- 
adjacent  subbands  is  comparable  to  stopband  attenuation.  Due  to  its  highly  symmetric  structure, 
the  PQF  can  be  implemented  much  more  efficiently  than  the  ordinary  wavelet  packet  transform 

m. 

The  difficulty  of  designing  a  pseudo-QMF  bank  lies  in  creating  a  single  low -pass  FIR  prototype 
filter  from  which  the  rest  of  the  filter  bank  can  be  constructed  [7]  [23]  [24].  Once  the  prototype  is 
obtained,  its  impulse  response  coefficients,  h(n),  are  modulated  by  a  cosine  term  to  created  each 
subband  filter.  The  f  bandpass  filter,  ,  of  an  M  channel  analysis  filter  bank  is  obtained  by 
[8]; 


H,in)  =  hin)cos(^n-i)  (2.27a) 
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Similarly,  the  synthesis  filters,  G,^{n),\  <  A:  <  M,  are  obtained  by  [8]; 


G,{n)  =  h{n)cos(^n  +  f )  (2.27b) 

The  operation  of  performing  M  filter  eonvolutions  with  a  window  size  of  N  requires  Mx  N 
multiplications  and  Mx(N-\^  additions.  The  implementation  of  the  PQF  requires  2M^  +  N 

multiplications  and  2M^  -  3M  +  N  additions  when  the  number  of  subbands  is  a  power  of  2  and 
the  prototype  filter  has  length  of  16  times  the  number  of  subbands  [8].  Figure  2.7  demonstrates 
the  computational  savings  of  the  pseudo-QMF  filter  bank,  when  using  the  previously  mentioned 
optimization  [8],  versus  performing  M  separate  filter  convolutions. 


Number  of  Multiplications  for  Subband  Decomposition  Number  of  Additions  for  Subband  Decomposition 


Number  of  Subbands  Number  of  Subbands 


Figure  2.7  Analysis  filter  bank  (a)  multiplications  (b)  additions 
Table  2.1  Analysis  filterbank  (a)  multiplications  (b)  additions 


4 

8 

16 

32 


The  filter  bank  is  built  from  a  single  prototype  filter  which  is  designed  to  satisfy  near-perfect 
reconstruction  conditions  within  a  certain  tolerance  [7].  The  design  of  such  a  prototype  filter 
generally  requires  solving  a  nonlinear  constrained  optimization  problem  [24].  For  our 
simulations,  we  used  the  prototype  filter  coefficients  specified  in  the  MPEG  audio  standard  [25]. 


Analysis  Filterbank  Additions 

subbands  s 

tandard 

PQF 

252  84 

1016  232 

3840  720 

15872  246 

1 

Analysis  Filterbank  Multiplications 

subbands  s 

tandard 

PQF 

4 

256  96 

8 

1024  256 

16 

4096  768 

32 

16384  2560 

2.7  Spectral  Characteristics  of  Digital  Linearly  Modulated  Signals 
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In  order  to  develop  an  efficient  bit  allocation  scheme  for  subband  coding,  it  is  important  to 
determine  the  spectral  content  of  the  digitally  modulated  signal.  Consider  a  random  information 
sequence  }  which  is  to  be  transmitted  over  a  channel.  A  modulator  maps  the  digital 
information  into  analog  waveforms  that  match  the  characteristics  of  the  channel  [10].  We  will 
consider  the  class  of  memoryless  linear  modulation  methods  with  the  following  base-band 
representation  [10]: 


00 

{12%) 

«=-oo 

Where  }  is  a  set  of  digital  symbols,  gj.(0  is  a  pulse-shaping  fdter  (Fig  2.8),  and  the  pulse 
duration  is  T. 


X 

Y^At-nT) 


v(0 


Figure  2.8  Digital  sequence  shaped  by  transmission  filter 


Examples  of  modulation  methods  belong  to  this  class  include  pulse  amplitude  modulation, 
phase-shift  keying,  and  quadrature  amplitude  modulation  [9]  [10].  Assuming  that  {/„}is  a  wide- 
sense  stationary  process  with  mean,  and  variance,  erf,  the  power  spectrum  density  (PSD)  of 
the  modulated  signal  (in  baseband  form)  can  be  expressed  as  [10]: 


CO 

=  +  A  YfiftAf-f)  (2.29) 

w=-oo 


I  i2 

Where  \G{f)\  is  the  PSD  of  the  pulse  shaping  filter,  gj{t) .  When  the  information  symbols  {/„ } 
are  equally  likely  and  symmetrically  distributed  in  the  complex  plane  [10],  //,  =  0,  and  the 
discrete  part  of  equation  2.29  vanishes: 

=  (2.30) 


This  implies  that  under  the  given  conditions,  the  PSD  of  the  baseband  signal  is  simply  a  scaled 

I  |2 

version  of  the  PSD  of  the  pulse  shaping  fdter,  \G{f)\  . 


It  is  assumed  that  the  sensors  we  will  be  using  are  provided  with  the  carrier  frequency  of  the 
signal,  so  that  we  can  remove  the  carrier  component  and  encode  the  complex  envelope  of  the 
signal.  Given  the  complex  envelope  of  a  digitally  modulated  signal  of  the  form  in  equation  2.28, 
assuming  that  the  sequence  {f }  is  a  zero-mean  stationary  random  process,  if  we  capture  a 
sufficiently  long  portion  of  the  signal  and  observe  its  spectrum,  we  can  compute  the  Fisher 
Information  and  determine  an  appropriate  bit  allocation.  The  bit  allocation  will  be  effective  for 
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the  remaining  duration  of  the  signal  so  long  as  the  pulse  shaping  filter  or  symbol  rate,  Yj,  does 
not  change. 
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Chapter  3 


The  Effects  of  Decimation 

3.1  Motivation  for  decimation 

In  order  to  accurately  perform  the  TDOA  estimation  process,  it  is  necessary  to  sample  the 
observed  continuous-time  signal  at  a  very  high  rate.  The  cross-correlation  of  two  discrete-time 
signals  will  yield  a  discrete  sequence;  therefore,  when  searching  for  the  peak  of  the  cross- 
correlated  signals,  the  precision  of  the  estimate  will  be  limited  to  integer  multiples  of  the 
sampling  period  (assuming  interpolation  is  not  used). 

Consider  a  signal  which  is  sampled  at  a  rate  of  20  MHz,  providing  a  resolution  of  50  ns.  If  we 
quantize  the  discrete-time  signal  using  only  one  bit  per  sample,  the  overall  rate  will  be  20  Mbps, 
which  is  infeasible  for  many  real-time  systems.  Additionally,  if  we  oversample  band-limited 
communication  signals,  it  is  possible  to  resample  at  the  Nyquist  rate  without  introducing 
significant  distortion.  In  an  effort  to  further  reduce  the  overall  code  rate,  we  will  explore  the  use 
of  decimation  to  reduce  the  sampling  rate  used  for  transmission. 


3.2  Overview  of  decimation 

When  the  continuous-time,  baseband  signal  is  sampled  at  a  rate  of  =  ^  ,  its  resulting 
spectrum  will  be  periodic  with  a  period  of  F^.  We  will  consider  the  region  centered  at  the  origin, 
confined  to  the  interval,  /  e  or  equivalently,  co  e  (Fig.  3.1). 


-.T  /T 

Figure  3.1  Spectra  of  a  discrete -time  signal 


The  down-sampling  of  a  signal  will  effectively  “expand”  [25]  its  spectra.  If  the  rate  is  reduced 
below  the  Nyquist  rate,  it  may  be  necessary  to  pass  the  signal  through  a  low-pass  filter,  H^{z), 
prior  to  down-sampling  in  order  to  avoid  aliasing  of  the  folded  spectrum,  (Fig.  3.2). 
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Figure  3.2  The  process  of  decimation 


At  the  receiver  end,  we  wish  to  restore  the  signal  to  its  original  sampling  rate  before  cross¬ 
correlation  in  order  to  regain  fine  precision  in  the  TDOA  estimate.  Up-sampling  the  signal  will 
“compress”  [25]  its  spectra,  causing  images  of  the  signal  to  be  introduced  into  the  upper  end  of 
the  spectrum.  These  images  can  be  removed  by  passing  the  up-sampled  signal  through  an 
interpolation  filter,  H^{z),  (Fig.  3.3). 


Figure  3.3  The  process  of  interpolation 


A  decimation/interpolation  codec  which  was  implemented  for  the  geo-location  scenario  is 
depicted  in  figure  3.4.  The  discrete-time  baseband  signal,  x{n),  is  first  passed  through  a  low- 
order  FIR  filter,  H^{z),  and  is  then  downsampled  by  a  factor  of  L.  When  the  decimation  factor, 
L,  is  very  large,  the  decimation  process  is  divided  into  two  stages.  The  decimated  signal  is  then 
transmitted  over  a  channel  to  a  receiver  which  up-samples  the  signal  by  the  factor,  L,  and  passes 
it  through  an  interpolation  filter,  H^{z).  The  reconstructed  signal,  x{n),  is  then  used  to  perform 
the  TDOA  estimation. 


Decimation  Interpolation 


— ► 

. > 

— ► 

Figure  3.4  Structure  of  a  decimation/interpolation  codec 


3.3  Experimental  Results 

In  our  experiments,  we  used  two  different  pulse-shaping  filters;  a  rectangle  (Fig.  3.5a),  and  a 
raised-cosine  (Fig.  3.5b). 
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time  time 
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Figure  3.5  (a)  rectangular  pulse  (b)  raised-cosine  pulse 

We  will  first  derive  the  power  spectrum  density  of  the  linearly  modulated  signals  using  the  two 
filters.  The  rectangular  pulse  is  defined  as; 

g(0  =  40<t<T  (3.1) 


Which  has  the  power  spectrum  density  (PSD); 


\G{ff  =  {AT)^  ■ 
V 


sin;^ 

Afr  J 


(3.2) 


Using  the  result  obtained  in  equation  2.30,  we  can  express  the  PSD  of  the  linearly  modulated 
baseband  signal  as; 


smT^i 

Afr 


(3.3) 


Let,  ^(/),  denote  the  PSD  of  the  band-limited  random  process,  sampled  at  the  rate  F^  =  ^  .  The 
PSD  of  the  sampled  signal  is  then; 


1  ” 


(3.4) 


s  «=-oo 


For  our  analysis,  we  will  consider  the  spectrum  on  the  interval  /  e 


</>{/)  =  ^CT^rT\ 


^  sin  Tfr 
V  Afr  y 


\f\^- 


(3.5) 
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Figure  3.6  compares  the  theoretically  obtained  PSD  (Eq.  3.5)  with  the  PSD  of  an  observed  signal 
where  the  pulse  duration  is  5ns,  and  it  is  sampled  at  a  rate  of  20MHz.  The  symbols,  {!„},  are 
i.i.d.,  uniformly  distributed  with  zero-mean. 


Normalized  frequency  magnitude 


frequency  / 

Figure  3.6  PSD  of  rectangular  pulse-shaped  signal  with  5  ns  period,  sampled  at  20  MHz 

The  PSD  of  the  rectangular  pulse  decays  proportional  to  as  indicated  in  equation  3.5.  Upon 
observation  of  the  PSD  (Fig.  3.6),  it  is  apparent  that  its  magnitude  is  very  small  in  the  range 
|/|  >  0.05,  which  means  that  the  signal  can  be  significantly  decimated  without  removing  a  large 
amount  of  energy.  In  order  to  better  illustrate  the  cumulative  distribution  of  power  in  the 
spectrum  of  a  continuous-time  rectangular  pulse,  we  have  numerically  integrated  its  PSD  (Fig. 
3.7). 

Cumulative  Power  Distribution 


Figure  3.7  Cumulative  power  distribution  of  a  rectangular  pulse 
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Notice  that  the  eumulative  power  grows  at  a  large  rate  for  low  frequeneies  (/  <  y),  and  at  high 
frequeneies,  the  eumulative  power  approaehes  1  asymptotically. 

To  illustrate  the  effects  of  decimation  on  the  spectrum  of  the  signal,  we  deeimate  an  observed 
signal  by  a  factor  of  50,  effectively  removing  all  but  the  central  lobe  (Fig.  3.8). 


Normalized  frequency  magnitude 
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frequency  / 

Figure  3.8  PSD  of  deeimated  signal 

After  the  signal  has  been  transmitted  to  a  fusion  center,  prior  to  cross-correlation,  the  signal  is 
upsampled  to  its  original  sampling  rate  and  is  then  interpolated.  Figure  3.9  shows  the  speetrum 
of  the  interpolated  signal  which  is  used  for  TDOA  estimation.  Notiee  that  the  eentral  lobe  is 
slightly  attenuated  at  its  edges  due  to  the  non-ideal  transition  bands  of  the  deeimation  and 
interpolation  filters. 


Normalized  frequency  magnitude 


frequency  /  F 


Figure  3.9  PSD  of  interpolated  signal 

The  second  pulse-shaping  filter  we  experimented  with  is  the  raised-cosine,  of  the  form; 
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8(0  =  y 


1 

1  +  COS  — 
T 


(  j\ 

t - 

,  2J 


0<t<T  (3.6) 


Which  has  a  power  spectrum  density  (PSD)  defined  as: 


\G{ft  = 


AT 


sin^  7^ 


{7fr)\\-fT^f 


(3.7) 


Using  equation  2.30,  we  ean  express  the  PSD  of  the  linearly  modulated  baseband  signal  as: 

cr^A^  sin^  Tfr 


(/>{/)=  -y  o  ,  , 

4  {7f)\i-fTy 

The  spectrum  of  the  sampled  signal,  eentered  at  /  =  0  is: 

ajA^  sin^  7^ 


(3.8) 


</>(/)  = 


47/  (T^Yii-rr) 


f\  <  W 
2rr2\2  AJ  I  “  /2 


(3.9) 


Figure  3.10  eompares  the  result  from  equation  3.9  to  the  PSD  of  an  observed  signal,  shaped  by  a 
raised-cosine  pulse  with  a  duration  of  5  ns,  sampled  at  20  MFlz.  The  modulated  symbols,  {/„ }, 
are  i.i.d.,  uniformly  distributed,  with  zero-mean. 


Normalized  frequency  magnitude 


frequency  /  F 


0.25 


Figure  3.10  PSD  of  raised  eosine  pulse  with  5  ns  period,  sampled  at  20  MHz 


Notice  that  the  central  lobe  of  the  speetrum  of  the  raised  cosine  is  wider  than  central  lobe  of  the 
sine  function  (it  does  not  have  zeros  at  /  =  ■)4’5//)-  Also,  the  PSD  of  the  raised  cosine  decays 
proportional  to  .  As  figure  3.10  demonstrates,  the  magnitude  of  the  side  lobes  is  very  small 

(not  even  visible  in  the  graph).  Figure  3.11  eontains  a  plot  of  the  eumulative  power  distribution 
of  a  raise-cosine  pulse,  similar  to  figure  3.7  for  the  rectangular  pulse.  Notice  that  the  eumulative 
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power  grows  at  a  large  rate  for  low  frequencies  (/  <j),  and  at  high  frequencies,  the  cumulative 
power  approaches  1  asymptotically. 


Cumulative  Power  Distribution 


Figure  3.11  Cumulative  power  distribution  of  a  rectangular  pulse 


To  illustrate  the  effects  of  decimation  on  the  spectrum  of  a  raised-cosine  shaped  signal,  we 
decimate  an  observed  signal  by  a  factor  of  25,  discarding  all  side  lobes  (Fig.  3.12). 


Normalized  frequency  magnitude 


Figure  3.12  PSD  of  decimated  signal 


Figure  3.12  shows  the  PSD  of  the  interpolated  signal  used  for  TDOA  estimation.  Notice  that  the 
PSD  of  the  interpolated  signal  is  not  very  different  than  the  original  signal,  as  the  majority  of  the 
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power  is  contained  in  the  central  lobe.  The  central  lobe  is  slightly  attenuated  at  its  edges  due  to 
the  non-ideal  transition  bands  of  the  decimation  and  interpolation  filters. 


Normalized  frequency  magnitude 


frequency  / 

Figure  3.13  PSD  of  interpolated  signal 

We  performed  a  series  of  experiments  to  investigate  the  effects  of  decimation  on  the  TDOA 
estimation  process.  In  the  experiments,  we  used  a  fixed  sampling  rate  of  20  MHz  and  generated 
signals  with  symbol  rates  up  to  2.5  million  symbols  per  second.  The  signals  were  comprised  of  a 
sequence  of  i.i.d.,  uniformly  distributed  symbols,  randomly  generated  with  zero-mean  unit 
variance,  and  modulated  using  a  pulse  with  duration,  T. 

The  sequence  length  and  symbol  rate,  were  jointly  chosen  to  produce  signals  with  a  duration 
on  the  order  of  3  ms.  To  simulate  centralized  estimation  (Fig.  1.2b),  where  two  sensors  share 
their  signals  with  a  fusion  center,  we  decimated  two  signals  by  a  factor,  L,  ranging  from  2  to  48, 
and  then  interpolated  by  the  same  factor  prior  to  cross-correlation.  The  resulting  TDOA  estimate 
was  compared  to  an  estimate  obtained  by  cross-correlating  the  same  signals  without  using 
decimation/interpolation.  Similarly,  to  simulate  distributed  estimation  (Fig  1.2a),  we  decimated 
and  interpolated  only  one  of  the  two  signals. 

Figures  3.14  and  3.15  show  the  variance  of  the  TDOA  estimation  error  obtained  from  performing 
the  decimation  experiment  100  times  with  a  rectangular  pulse-shaping  filter.  The  results  in 
figure  3.14  are  for  the  centralized  estimation  case,  figure  3.15  depicts  the  results  of  the 
distributed  estimation  case.  In  the  graphs,  the  symbol  rate,  relative  to  the  sampling  rate  is  given 
in  the  units  samples  per  symbol. 


23 


TDOA  6rror  g  TDOA  error 


0  400 

Decimation  Factor  Samples  per  Symbol 

Figure  3.14  TDOA  error  as  a  function  of  symbol  rate  and  decimation  factor 

(rectangle,  centralized) 


Decimation  Factor  Samples  p^^  gymbol 

Figure  3.15  TDOA  error  as  a  function  of  symbol  rate  and  decimation  factor 

(rectangle,  distributed) 
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In  comparing  figure  3.14  with  figure  3.15,  it  is  apparent  that  the  TDOA  estimation  error  is 
significantly  less  when  only  one  of  the  two  signals  is  decimated  and  interpolated.  Additionally, 
in  both  cases,  the  TDOA  error  introduced  by  decimating  signals  with  low  symbols  rates  is 
minimal. 


The  same  decimation  experiments  were  also  performed  with  signals  which  were  shaped  using  a 
raised-cosine  pulse.  Figures  3.16  and  3.17  depict  the  results  of  the  raised-cosine  experiments. 
Once  again,  the  TDOA  error  introduced  by  decimation  in  the  distributed  estimation  scenario 
(Fig.  3.17)  is  significantly  less  than  the  error  introduced  in  the  centralized  estimation  scenario 
(Fig.  3.16). 
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Figure  3.16  TDOA  error  as  a  function  of  symbol  rate  and  decimation  factor 

(raised-cosine,  centralized) 
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Figure  3.17  TDOA  error  as  a  function  of  symbol  rate  and  decimation  factor 

(raised-cosine,  distributed) 


Comparing  the  results  of  the  decimation  experiments  for  the  rectangular  pulse  (Fig.  3.14)  with 
those  for  the  raised-cosine  pulse  (Fig.  3.16),  it  can  be  seen  that  for  low  decimation  factors,  the 
raised-cosine  shaped  signals  produce  less  TDOA  error  than  the  rectangular  shaped  signals. 
However,  for  large  decimation  factors,  the  TDOA  error  for  the  rectangular  shaped  signals  is 
slightly  less  than  for  the  raised-cosine  shaped  signals.  Figures  3.18a  and  3.18b  compare  the 
TDOA  error  variance  for  both  cases  when  the  signals  are  decimated  by  the  factors  10,  20,  30, 
and  42. 

Rectangle  Raised  Cosine 


Figure  3.18  TDOA  error  as  a  function  of  symbol  rate  for  decimation  factors  10,  20,  30,  42 
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(a)  rectangular  shaped  signal  (b)  raised-cosine  shaped  signal 

The  reason  for  this  behavior  may  be  explained  by  examining  the  spectral  shape  of  the  pulses. 

The  PSD  of  the  raised-cosine  decays  proportional  to  therefore,  the  power  contributed  by 

side  lobes  is  less  than  the  power  contributed  by  the  side  lobes  of  the  PSD  of  the  rectangle  (which 
decays  proportional  to  Also,  the  central  lobe  of  the  PSD  of  the  raised-cosine  is  wider  than 

the  central  lobe  of  the  PSD  of  the  rectangle.  Figure  3.19  depicts  the  amount  power  which  is 
retained  after  decimating  a  rectangular  pulse  and  a  raised-cosine  pulse. 


Signal  Power  as  a  Function  of  Decimation  Factor 


Figure  3.19  Signal  power  as  a  function  of  decimation  factor 
(Decimation  removes  energy  from  the  signal) 


When  we  decimate  by  small  factors  (but  large  enough  to  reduce  the  rate  below  the  Nyquist  rate), 
we  are,  in  effect,  removing  side  lobes  from  the  signal;  therefore,  more  energy  will  be  removed 
from  the  signal  shaped  by  a  rectangular  pulse.  As  a  result,  the  signal-to-noise  ratio  of  the 
interpolated  rectangular  pulse-shaped  signal  will  be  lower,  which  will,  in  turn,  have  a  greater 
effect  on  the  estimation  error. 

When  we  decimate  by  large  factors  (large  enough  to  remove  all  side  lobes),  we  remove  power 
from  the  central  lobe.  The  majority  of  the  power  in  a  raised-cosine  pulse  lies  in  its  central  lobe, 
in  the  region  |/|  <  y ,  while  the  majority  of  the  power  in  the  rectangular  pulse  is  in  the  central 
lobe,  in  the  region  |/|  <  y .  Decimating  the  sampled,  raised-cosine  shaped  signal  by  a  factor, 
L>y^,  will  result  in  removing  power  from  the  central  lobe;  whereas,  the  central  lobe  of  the 
rectangular  shaped  signal  will  not  be  affected  until  we  decimate  by  a  larger  factor,  T  >  ^ . 
Notice  (Fig.  3.19)  that  the  rectangular  shaped  signal  is  decimated  by  factors  greater  than  about 
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,  it  retains  more  power  than  the  raised-eosine  shaped  signal,  decimated  by  the  same  amount. 
This  loss  of  power  at  higher  decimation  factors  is  what  accounts  for  the  improvement  of  the 
TDOA  estimate  for  the  rectangular  shaped  signal  over  the  raised-cosine  shaped  signal. 
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Chapter  4 


The  Effects  of  Quantization 

4.1  The  Subband  Encoder 

In  this  chapter  we  will  follow  the  work  of  [19-22]  in  developing  a  subband  encoder  which  is 
designed  to  compress  an  intercepted  signal  while  minimizing  adverse  effects  on  the  TDOA 
estimation  process.  The  observed  signal  is  decomposed  into  a  number  of  frequency  subbands, 
each  of  which  is  quantized  independently,  using  a  uniform  scalar  quantizer.  The  structure  of  the 
subband  codec  is  depicted  in  figure  4. 1 .  The  sequence  that  is  being  compressed,  x{n),  is  the 
sampled  complex  envelope  of  an  observed  signal,  which  is  assumed  to  be  finely  quantized  using 
32  bits  per  sample.  The  real  and  imaginary  parts  of  x{n)  are  independently  quantized;  therefore, 
the  compression  ratio  of  the  subband  encoder  can  be  calculated  as  CR  =  ^327? ,  where  R  is  the 
number  of  bits  budgeted  per  sample  for  the  compressed  signal. 
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Figure  4.1  Subband  encoder/decoder  block  diagram 


*-x{n) 


The  subband  encoder  consists  of  anM-channel  Pseudo-QMF  analysis  filter  bank  (Section  2.6) 
and  M  uniform  scalar  quantizers.  The  bit  allocation  for  the  quantizers  is  determined  by  a  fusion 
center  which  observes  a  portion  of  an  intercepted  signal.  The  focus  of  this  chapter  will  be  on 
comparing  the  effectiveness  of  a  bit  allocation  determined  using  Fisher  Information-based 
criteria  versus  choosing  an  allocation  based  solely  on  mean-squared-error  distortion  criteria. 


4.2  Mean  Squared  Error  Distortion  Criteria 

In  order  to  establish  a  basis  for  comparison,  we  will  first  use  the  well-known  MSB  distortion 
criteria  to  determine  a  bit  allocation  for  the  subband  encoder.  We  use  the  same  rate-distortion 
model  as  in  [19]  for  the  subband; 

A(^,)  =  ,\<i<M  (4.1) 
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where  the  variance  of  the  signal  samples  of  the  i‘^  subband,  erf,  is  determined  empirically  from  a 
portion  of  the  signal.  Using  the  M  rate-distortion  functions  and  a  rate  constraint,  R,  the  bit 
allocation,  B,  is  determined  by  solving  the  constrained  integer  optimization  problem; 

,<i?|  (4.2) 

Where  the  bit  allocation  vector  is  B  =  \b.  >  0|z  e  1,2,-  • 

To  demonstrate  the  bit  allocation  process,  we  will  use  a  rectangular  pulse-shaped  signal  (Section 
3.3)  with  a  pulse  duration,  T  =  5ns,  sampled  at  the  rate,  =  20 MHz .  Since  the  band-limited 

communication  signal  is  oversampled,  it  can  be  downsampled  prior  to  subband  encoding  without 
introducing  significant  error  in  the  TDOA  estimate  (Chapter  3).  The  signal  is  downsampled  by  a 
factor  of  16  so  that  the  central  lobe  and  two  side-lobes  of  the  power  spectrum  density  are 
retained.  The  magnitude  spectrum  of  the  downsampled  signal  is  shown  in  figure  4.2. 


t  to; 
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Normalized  frequency  magnitude 


Figure  4.2  Magnitude  spectrum  of  a  rectangular  shaped  signal 
(T=5ns,  Fs=20MHz,  SNR=10dB) 


The  signal  is  decomposed  into  32  frequency  subbands  using  a  pseudo-QMF  bank.  Figure  4.3 
depicts  the  variance  of  the  subband  samples,  cr,^,l  <  z  <  32. 
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Variance  of  Subband  Samples 
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Figure  4.3  Variance  of  the  subband  samples 


2 

The  variances,  cr,  ,1  <  z  <  32,  are  then  used  in  equation  4.1  to  compute  the  rate-distortion 
functions  for  the  32  subbands.  The  rate-distortion  functions,  combined  with  the  rate  constraint, 
R,  allow  us  to  formulate  the  constrained  optimization  problem  in  equation  4.2.  The  optimization 
problem  was  solved  using  a  MATLAB  implementation  of  the  Lagrange  optimization  algorithm 
described  in  Appendix  A.  Figure  4.4  depicts  the  32  rate-distortion  functions  and  the  solution 
points,  b.,l  <  i  <  M,  (redx’s)  when  the  rate  is  constrained  to  4  bits  per  sample  (4:1  compression 
ratio). 


Figure  4.4  Rate-Distortion  functions  for  the  32  quantizers  (4:1  Compression  Ratio) 

The  bit  allocation,  B,  which  provides  a  solution  to  the  constrained  optimization  problem  (Eq. 
4.2),  for  this  example  is  depicted  in  figure  4.5 
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Subband  Bit  Allocation 
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Figure  4.5  Bit  allocation  computed  using  MSB  criteria 

Notice  that  by  using  the  MSB-based  distortion  criteria,  bits  are  allocated  to  subbands  based 
primarily  on  their  power  content.  Since  white  noise  was  added  to  the  signal  prior  to 
compression,  the  PSD  of  the  noisy  signal  at  high  frequencies  is  nearly  flat;  therefore,  the  bits  are 
allocated  almost  evenly  over  the  side  lobes. 


4.3  Fisher  Information  Distortion  Criteria 


The  Bisher  Information-based  distortion  that  we  will  use  is  that  developed  in  [19-22]  (Section 
2.5).  We  will  derive  the  Bisher  Information  (for  TDOA)  for  the  rectangular  and  raised  cosine 
pulses  (section  3.3)  using  the  formulas  introduced  in  section  2.6.  Continuing  with  the  power 
spectrum  density  of  a  rectangular  pulse-shaped  signal,  sampled  at  the  rate  (Eq-  3.5); 
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\f\^- 


(4.3) 


Substituting  the  PSD  of  the  signal  into  the  Bisher  Information  formula  that  was  presented  earlier 
(Bq.  2.20),  we  have; 
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(4.4) 


The  terms  in  the  numerator  and  denominator  of  the  integrand  of  equation  4.4  cancel  and  the 
Bisher  Information  becomes; 
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(4.6) 
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As  discussed  in  section  2.6,  in  the  wireless  sensor  network  scenario,  it  is  assumed  that  we  cannot 
observe  the  original  signal,  uncorrupted  by  additive  noise.  Additionally,  it  is  assumed  that  the 
sensor  which  is  performing  the  compression  will  know  the  variance  of  the  noise  at  a  neighboring 
sensor;  therefore,  it  is  not  feasible  to  directly  compute  the  Fisher  Information  (Eq.  4.6).  Instead, 
we  follow  the  approach  of  [19]  [22]  and  approximate  the  Fisher  Information  using  the  DFT 
coefficients  of  an  intercepted,  noisy  signal  (Eq.  2.25).  To  demonstrate  the  performance  of  the 
approximate  FI  model  (Eq.  2.25),  figure  4.6  compares  the  quadratically  weighted  DFT 
coefficients  of  an  intercepted  noisy  signal  to  the  integrand  of  equation  4.6  (where  the  integrand  is 
sampled  using  the  same  number  of  points  as  the  DFT). 


Quadratically  Weighted  Spectrum 


Figure  4.6  Quadratically  weighted  power  spectrum  of  rectangular  pulse  shaped  signal 


Because  the  FI  is  proportional  to  the  integral  of  the  quadratically  weighted  PSD,  we  can  examine 
the  contribution  of  individual  frequency  components  to  the  Fisher  Information.  Due  to  the  form 
of  the  FI  of  the  rectangular  pulse-shaped  signal  (Eq.  4.6),  it  is  apparent  that  the  side  lobes 
contribute  significantly  to  the  FI. 


Similarly,  for  the  raised-cosine  pulse  shaped  signal  (Eq.  3.6),  it  can  be  shown  that  the  Fisher 
Information  with  respect  to  TDOA  estimation  is; 
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df 


(4.7) 


Figure  4.7  depicts  the  integrand  of  equation  4.7  and  the  quadratically  weighted  DFT  coefficients 
of  an  intercepted,  noisy  raised-cosine  pulse-shaped  signal.  Notice  that  the  weighted  spectrum 
decays  proportional  to  ^4 ,  implying  that  frequency  components  in  the  range  |/|  <  y  will 

contribute  heavily  to  the  Fisher  Information. 
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Quadratically  Weighted  Spectrum 


Following  the  work  of  [19-22],  we  use  the  Fisher  Information-based  model  (Eq.  2.25,  2.26)  to 
characterize  the  effects  of  quantization  on  the  TDOA  estimation  process.  The  FI  corresponding 
to  the  subband  for  the  multi-bit  case  is  computed  as  [19] [21]; 


Ncr^  +  q] 


(4.8) 


The  distortion  introduced  by  quantization  of  the  subband,  q],  is  calculated  using  the  rate- 
distortion  model  in  equation  4. 1 .  The  FI  corresponding  to  the  subband  for  the  single-bit 
quantization  case  is  computed  as  [19]: 
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Where  A,  is  the  set  of  indices  belonging  to  the  subband.  The  bit  allocation  is  then  determined 
by  solving  the  optimization  problem; 

(4.10) 


subject  to; 


M 
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Where  the  bit  allocation  vector,  B  =  >  0|z  e  1,2,-  •  -,4/ j,  is  constrained  to  non-negative  integers 

and  the  objective  function  is: 
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0A=0 
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(4.11) 


A  Lagrange  optimization  routine  was  implemented  in  MATLAB  to  solve  the  eonstrained 
minimization  problem  presented  in  seetion  4.2  (Eq.  4.2).  If  we  eould  somehow  reformulate  the 
eonstrained  maximization  problem  (Eq.  4.10)  as  a  minimization  problem,  then  we  eould  use  the 
same  optimization  proeedure  to  solve  for  both  the  MSE-based  and  El-based  bit  alloeations.  The 
objeetive  funetion,  J.,  is  defined  by  3  different  funetions  (Eq.  4.1 1),  depending  on  the  value  of 
the  rate,  b..  Consider,  first,  the  funetion,  j\”'\  for  multi-bit  quantization  (Eq.  4.8).  Replaeing  the 
numerator  of  equation  4.8  by  G,  ean  be  rewritten  as: 
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Using  equation  4.12b,  the  maximization  problem  (Eq.  4.10),  for  the  ease  when  b.  >  2  ean  be 
expressed  as: 
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Where  B  -  \b.  >  2\i  e  1,2,-  •  Sinee  q],  by  definition  (Eq.  4.1),  is  the  only  term  in  equation 
4.13  whieh  is  a  funetion  of  b^,  the  parameter  over  whieh  we  are  searehing,  the  optimization 
problem  (Eq.  4.13)  ean  be  equivalently  stated  as  [22]: 
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Where  B  =  {h,  >  2\i  e  1,2,-  •  Similarly,  we  reformulate  the  values  of  the  objeetive  funetion 
of  the  maximization  problem  (Eq.  4.1 1)  for  the  single-bit  (h,  =  1)  and  zero-bit  (h,  =  0)  oases: 
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The  bit  allocation  can  now  be  determined  by  solving  the  optimization  problem; 


min' 

B 


subject  to: 


Where  B  =  {Z),  >  0|/  e  1,2,-  •  •,M^  and  the  objective  function  is  defined  as; 
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Figure  4.8  depicts  the  objective  functions  (Eq.  4.17)  computed  for  a  32  subband  encoder.  The 
red  x’s  denote  the  solution  points,  b.,\  <  z  <  M,  to  the  optimization  problem  when  the  rate  was 
constrained  to  4  bits  per  sample  (4:1  compression  ratio). 


X  10* 


Figure  4.8  FI-based  objective  function  for  the  32  quantizers  (4:1  Compression  Ratio) 

Notice  that  the  objective  functions  in  figure  4.8  are  not  convex  due  to  the  values  given  by  the 
single-bit  approximation.  As  a  result  of  the  non-convexity,  the  Fagrange  optimization  algorithm 
(section  2.4)  used  to  solve  the  constrained  optimization  problem  (Eq.  4.17)  may  not  yield  an 
optimal  solution.  Instead,  we  use  the  Fagrange  optimization  procedure  to  find  the  closest 
solution  from  below  (a  rate  which  is  less  than  or  equal  to  the  constraint),  then  we  use  a  water¬ 
filling  method  (Appendix  A)  to  allocate  any  remaining  bits.  Figure  4.9  depicts  the  bit  allocation 
determined  by  solving  the  constrained  optimization  problem  (Eq.  4.17). 
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Subband  Bit  Allocation 


1  4  8  12  16  20  24  28  32 

Subband  Index 


Figure  4.9  Bit  allocation  computed  using  FI-based  criteria 

Notice  that  more  bits  are  allocated  to  the  side  lobes  when  we  use  the  FI-based  approach  (due  to 
quadratic  weighting  of  the  spectrum).  In  the  MSE-based  approach,  the  bits  were  allocated  more 
evenly  in  the  subbands  which  contain  the  side  lobes. 


4.4  Experimental  Results 

We  conducted  a  series  of  subband  encoding  experiments  to  compare  the  performance  of  using  a 
bit  allocation  determined  via  MSE-based  criteria  versus  a  bit  allocation  determined  by  the  EI- 
based  approach.  In  the  experiments,  linearly  modulated  signals  were  randomly  generated  and 
had  an  average  duration  of  3  ms  (60,000  samples  with  =  20 MHz).  The  signals  were 

decimated  prior  to  subband  encoding  to  enhance  the  TDOA  error  introduced  by  fine  quantization 
(as  a  coarsely  quantized  signal  can  still  provide  an  error-free  TDOA  estimate  when  the  number 
of  samples  is  very  large).  Figure  4.10  depicts  the  codec  that  was  used  in  the  experiments. 


Subband  Encoder/Decoder 

Figure  4.10  Codec  block  diagram 

We  selected  the  symbol  rates  and  decimation  factors  using  data  obtained  from  the  decimation 
experiments  (Chapter  3).  Figure  4.1 1  depicts  3  operating  points  that  were  chosen  for  the 
experiments  with  a  rectangular  pulse-shaped  signal. 
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Table  4.1  Operating  points  for  subband  encoding  experiments 


Operating 

Point 

Symbol  Rate 
(samples/symbol) 

Decimation  Factor 

1  100 

48 

2  19 

14 

3  9 

8 

Samples  per  Symbol 


Figure  4.1 1  Operating  points  for  subband  encoding  experiments 
(rectangular  pulse-shaped  signal) 


The  red  regions  in  figure  4.1 1  indicate  large  TDOA  error  whereas  the  dark  blue  regions  indicate 
very  small  TDOA  estimation  error.  The  operating  points  were  chosen  to  be  in  regions  where 
decimation  produced  little  to  no  error,  but  where  further  decimation  would  result  in  a  sizeable 
increase  in  error. 


In  the  first  set  of  experiments,  we  randomly  generated  100  signals  and  independently  computed  a 
bit  allocation  for  each  received  signal.  A  fixed  compression  ratio  was  chosen  and  the  SNR  of  two 
signals  (received  by  two  different  sensors)  was  varied  from  40dB  to  lOdB.  The  scenario  of 
centralized  estimation  (where  two  sensors  compress  their  observed  signals  prior  to  transmitting 
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them  to  a  fusion  center  for  processing)  was  simulated.  The  compressed/decompressed  signals 
were  used  to  obtain  a  TDOA  estimate,  which  was  then  compared  to  a  TDOA  estimate  obtained 
from  signals  which  had  not  been  compressed  (then  decompressed). 


Subband  Encoding  Petformance,  100  sam/sym,  CR=8;1 ,  F=48  Subband  Encoding  Performance.  100  sam/sym,  CR=4:1 ,  F=48 


Figure  4.12  Subband  encoder  performance  using  bit  allocation  for  each  signal 
(a)  exp  1,  CR=8:1  (b)  exp  1,  CR=4:1  (c)  exp  2,  CR=8:1  (d)  exp  2,  CR=4:1 
(e)  exp  3,  CR=8;1  (f)  exp  3,  CR=4;1 
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Figures  4.12  and  4.14  depict  the  results  of  the  experiment  performed  at  the  3  different  operating 
points  (table  4.1).  Figures  4.12  (a)  and  (b)  (top)  present  the  results  for  operating  point  1,  figures 
4.12  (c)  and  (d)  (middle)  present  the  results  for  operating  point  2,  and  figures  4.12  (e)  and  (f) 
(bottom)  present  the  results  for  operating  point  3.  Figures  4.12  (a),(c),  and  (e)  (left)  depict  the 
results  when  the  compression  ratio  was  fixed  at  8:1,  and  Figures  4.12  (b),(d),  and  (f)  (right) 
depict  the  results  when  the  compression  ratio  was  fixed  at  4: 1 . 

Notice  (Fig.  4.12)  that  subband  encoding  using  the  Fl-based  approach  resulted  in  less  TDOA 
estimation  error  than  the  MSE-based  approach.  According  to  the  approximate  FI  model  that  was 
used  for  compression  (Eq.  4.8,  4.9),  the  El  should  increase  as  the  signal-to-noise  ratio  of  the 
received  signal  is  increased,  meaning  that  the  TDOA  error  variance  should  decrease.  Observing 
the  experimental  results  (Eig.  4.12),  it  is  apparent  that,  in  general,  the  TDOA  estimation  error  is 
decreasing  as  the  SNR  is  increased  up  to  about  25dB.  In  most  cases  (Eig.  4.12(a-e)),  the 
estimation  error  increases  as  the  SNR  is  increased  above  25  dB.  Referring  back  to  the 
approximate  El  model  (Eq.  4.8,  4.9),  we  mentioned  that  a  necessary  condition  for  the  model  to 
hold  was  that  the  variance  of  the  noise  introduced  by  quantization  remains  smaller  than  the 
variance  of  the  sensor  noise.  When  the  signal-to-noise  ratio  becomes  large,  the  variance  of  the 
sensor  noise  becomes  small  and  the  approximate  El  model  fails  to  hold. 

The  approach  of  calculating  a  bit  allocation  for  each  received  signal  can  be  computationally 
expensive  for  wireless  sensors  and  also  requires  additional  side  information  to  be  included  with 
each  compressed  signal.  We  subsequently  explored  a  suboptimal  approach  which  would  not 
require  computation  of  a  bit  allocation  for  each  received  signal.  In  the  second  set  of 
experiments,  100  signals  were  generated  at  random,  as  before,  but  only  one  of  the  signals  was 
analyzed  to  determine  the  bit  allocation  which  would  be  used  to  compress  all  100  signals.  To 
demonstrate  the  selection  of  a  bit  allocation,  figure  4.13(a)  compares  the  average  of  the 
allocations  calculated  for  100  randomly  generated  signals,  using  El  criteria,  to  the  allocations 
calculated  for  two  randomly  selected  message  signals  4.13(b)  and  (c). 
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Rate  (bits)  Rate  (bits)  Rate  (bits) 


Subband  Index  (a) 

FI  bit  allocation  message  #5 


Subband  Index  (b) 

FI  bit  allocation  message  #41 


Subband  Index  (c) 


Figure  4.13  FI -based  bit  alloeation  (a)  averaged  over  100  messages 
(b)  message  #5  (e)  message  #41 

Figure  4.14  depiets  the  results  of  the  experiment  using  the  suboptimal  approaeh,  performed  at  the 
3  operating  points  (table  4.1). 
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Subband  Encoding  Performance.  100  sam/sym,  CR=8:1 ,  F=48 


Subband  Encoding  Performance,  100  sam/sym.  CR=4;1 ,  F=48 


Figure  4.14  Subband  encoder  performance  using  one-time  bit  allocation 
(a)  exp  1,  CR=8;1  (b)  exp  1,  CR=4;1  (c)  exp  2,  CR=8;1  (d)  exp  2,  CR=4:1 
(e)  exp  3,  CR=8;1  (f)  exp  3,  CR=4;1 


Figures  4.14  (a)  and  (b)  (top)  present  the  results  for  operating  point  1,  figures  4.14  (c)  and  (d) 
(middle)  present  the  results  for  operating  point  2,  and  figures  4.14  (e)  and  (f)  (bottom)  present 
the  results  for  operating  point  3.  Figures  4.14  (a),(c),  and  (e)  (left)  depict  the  results  when  the 
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compression  ratio  was  fixed  at  8:1,  and  Figures  4.14  (b),(d),  and  (f)  (right)  depiet  the  results 
when  the  eompression  ratio  was  fixed  at  4: 1 . 

It  ean  be  seen  (Fig.  4.14)  that,  in  general,  the  Fl-based  approaeh  results  in  less  TDOA  estimation 
error  than  the  MSE-based  approaeh.  In  most  oases,  the  error  deoreases  as  the  SNR  is  inoreased 
until  about  25dB.  As  the  SNR  is  inoreased  above  25dB,  the  error,  in  some  oases,  inoreases 
(similar  to  the  first  set  of  experiments).  Notioe  in  figure  4.14(b)  that  the  error  seems  to  fluotuate 
as  the  SNR  in  inoreased.  This  behavior  is  likely  a  results  of  using  the  suboptimal  approaeh; 
however,  also  note  that  the  error  is  nearly  zero. 

In  oomparing  the  results  from  the  first  set  of  experiments  (Fig.  4.12)  to  the  results  of  using  the 
suboptimal  approaeh  (Fig.  4.14),  it  ean  be  seen  that  in  some  oases  (b,  o,  and  d)  that  the 
suboptimal  approaeh  aotually  produoed  less  error.  For  the  other  oases  (a,  e,  f),  the  inorease  in 
error  resulting  from  the  suboptimal  approaeh  was  reasonably  small. 

Upon  observing  the  results  of  the  experiments  whioh  were  performed,  (Fig.  4.12,  4.14),  one  ean 
oonolude  that  the  approximate  FI  model  (Eq.  4.8,  4.9)  oharaoterized  the  effeots  of  quantization  on 
the  estimation  task  reasonably  well.  Additionally,  the  suboptimal  approaeh  of  using  a  fixed  bit 
allooation  to  oompress  an  ensemble  of  messages  produoed  good  results. 
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Chapter  5 

Conclusion 

5.1  Discussion 

In  this  thesis,  we  explored  the  effects  of  compression  on  the  task  of  estimating  time-difference- 
of- arrival  parameters.  We  began  by  investigating,  experimentally,  the  effects  of  decimation  on 
TDOA  estimation  for  digital  linearly  modulated  signals  with  different  symbol  rates,  using  both 
rectangular  and  raised-cosine  pulse-shaping  transmit  filters.  We  showed  that  the  band-limited 
communication  signals  that  were  used  could  be  decimated  without  introducing  considerable 
error,  so  long  as  the  power  within  the  central  lobe  of  the  power  spectrum  was  retained. 

We  next  explored  the  effects  of  quantization  on  the  estimation  task  using  a  Fisher  Information- 
based  model  developed  in  [18-22].  Due  to  the  structure  of  the  FI  model  (Eq.  4.8)  we 
implemented  a  subband  encoder  which  would  allow  us  to  selectively  quantize  frequency 
subbands  based  on  their  contribution  to  the  Fisher  Information.  We,  additionally,  developed  a 
suboptimal  algorithm  to  perform  subband  encoding  by  analyzing  only  a  portion  of  the  signals 
that  were  to  be  compressed.  Experimental  results  showed  that  subband  encoding  using  the  EI- 
based  approach  resulted  in  less  TDOA  estimation  error  than  using  the  standard  MSE-based 
distortion  criteria. 


5.2  Further  work 

In  order  to  characterize  the  effects  of  quantization  on  the  TDOA  estimation  process,  we  used  a 
Eisher  Information-based  model  (Eq.  4.8,  4.9)  developed  in  [18-22].  The  model  relies  on  the 
assumptions  that  1)  the  sensor  noise  is  i.i.d.  AWGN,  2)  the  noise  introduced  by  quantization  is 
uniformly  distributed,  white,  and  independent  of  the  sensor  noise  [19],  and  3)  the  variance  of  the 
sensor  noise  is  much  larger  than  the  variance  of  the  quantization  noise.  When  the  intercepted 
signal  is  coarsely  quantized  or  its  SNR  is  high,  conditions  2  and  3  will  be  violated  and  the  El 
model  (Eq.  4.8)  will  not  hold  [22].  In  [19],  the  authors  provide  an  approximate  model  for  the 
case  of  single-bit  quantization  (Eq.  4.9).  Under  high  SNR  conditions,  the  author  of  [22]  provides 
a  method  for  approximating  the  Eisher  Information  which  is  of  large  complexity,  and  is  likely 
not  feasible  for  use  in  resource  constrained,  real-time  tracking  scenarios.  There  remains  further 
research  in  developing  a  El  model  of  low  complexity  which  can  hold  under  high  SNR  conditions. 

Similar  to  the  El-based  model  developed  for  quantization  [18-22],  one  could  explore  models  to 
characterize  the  effects  of  decimation  and  subsequent  interpolation  on  the  TDOA  estimation 
process. 

Another  possible  avenue  of  research  would  be  to  compare  the  performance  of  decimation  to  the 
performance  of  quantization  for  compressing  signals  for  TDOA  estimation.  In  addition  to 
analyzing  the  effects  on  the  estimation  task,  one  could  include  an  analysis  of  the  computational 
complexity  of  both  methods. 
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Continuing  under  the  framework  of  subband  encoding  for  TDOA  estimation,  one  could  perform 
further  analysis  to  determine  the  importance  of  different  frequency  subbands  to  the  estimation 
task.  In  this  thesis,  we  used  an  FI-based  model  which  essentially  applied  a  quadratic  weighting 
to  the  spectrum  of  the  signal  which  was  to  be  compressed.  Perhaps,  with  further  analysis,  one 
could  determine  the  role  which  different  frequency  components  play  on  the  estimation  of  the 
TDOA  parameter. 
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Appendix  A:  Lagrange  optimization  algorithm 


Given  a  set  of  diserete  (monotonieally  deereasing)  rate-distortion  funetions,  D.{K),\  <i<  M ,  and 
a  rate  eonstraint  R,  the  Lagrange  optimization  algorithm  will  determine  an  optimal  way  to 
distribute  a  budgeted  set  of  bits  to  a  set  of  quantizers  .  The  optimization  problem  whieh  we  are 
trying  to  solve  is: 


iM  I  M 

XA(^,)r  subjee  tto:^^,^^  (A.l) 

where  B  =  {b,  >  0|z  e  1,2,-  •  The  set  of  admissible  bit  alloeation  values,  b,  ,  is  eonstrained  to 
non-negative  integers. 

Introdueing  a  Lagrange  multiplier,  /I  >  0, 
uneonstrained  problem: 

min' 

B 

is  also  the  solution  to  the  budget-eonstrained  problem  in  equation  A.  1  Moreover,  they  show  that 
for  a  given  X : 


eferenees  [4] [16]  show  that  the  solution  to  the 


XA(^,)+Ab,  (A.2) 


iM  I  M 

X  A  A,)  +^i\  =  Zmin{A(^,)  +  }  (A.3) 

1=1  J  1=1 

As  a  result,  the  minimum  of  the  “Lagrangian  eost  funetion”,  D.{b.)  +  Ab, ,  ean  be  eomputed 
independently  for  eaeh  quantizer.  Figure  A.l  presents  a  graphieal  interpretation  of  minimizing 
the  Lagrangian  eost  funetion. 
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Convex  Hull 


Figure  A.l  Graphical  interpretation  of  minimizing  the  Lagrangian  cost  function 

As  the  rate,h.,  increases,  the  distortion,  D.(b.),  decreases;  therefore,  we  have  a  trade-off  between 
rate  and  distortion.  The  Lagrange  multiplier  allows  us  to  select  specific  trade-off  points. 
Minimizing  the  Lagrangian  cost  function  when  /I  =  0  is  equivalent  to  minimizing  the  distortion; 
and  minimizing  the  Lagrangian  cost  when  /I  ^  oo  is  equivalent  to  minimizing  the  rate.  As  can 
be  seen  in  figure  A.  1 ,  the  point  in  the  rate-distortion  function  of  the  quantizer  which 
minimizes  the  Langrangian  cost  is  that  point  at  which  the  line  of  absolute  slope  A,  is  tangent  to 
the  convex  hull  of  the  rate-distortion  function. 

The  optimization  algorithm  iteratively  searches  for  a  Lagrange  multiplier.  A,  which  yields  a 
solution  that  meets  the  predetermined  rate  constraint.  The  choice  of  an  initial  Lagrange 
multiplier,  A. ,  will  greatly  affect  the  rate  of  convergence;  therefore,  A.  was  chosen  as  a  function 
of  the  rate  constraint,  R.  We  used  the  algorithm  discussed  in  [16]  to  search  for  a  value.  A,  which 
meets  the  rate  constraint.  Figure  A.2  is  a  flow  diagram  of  the  algorithm  which  was  implemented 
(using  the  notation  from  [16]).  When  an  exact  solution  cannot  be  found  using  the  Lagrange 
optimization  alone,  the  Lagrange  optimization  is  used  to  find  the  closest  solution  from  below, 
then  a  water-filling  method  is  used  to  allocate  any  remaining  bits  (Approximate  Solution  at 
bottom-right  of  Fig.  A.2). 
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Figure  A.2  Flowchart  of  Lagrange  optimization  algorithm  (with  water-filling) 


The  Lagrange  multiplier-based  optimization  algorithm  will  provide  an  optimal  solution  so  long 
as  the  discrete  rate-distortion  functions  are  convex  and  contain  points  which  provide  a  solution 
for  a  given  rate  constraint  (i.e.  if  the  set  of  admissible  bit  allocation  values,  h,  ,  is  further 
constrained,  an  optimal  solution  may  be  unreachable). 
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Appendix  B:  Additional  Formulations 

B.l  Cramer-Rao  bound  for  time  delay  estimation 

This  appendix  gives  an  outline  of  the  derivation  of  the  Cramer-Rao  bound  for  time  delay 
estimation.  For  a  detailed  formulation,  the  reader  is  direeted  to  the  sourees  from  which  the 
outline  was  created  [4] [15] [16] [17]. 

The  Cramer-Rao  inequality  expresses  an  asymptotic  lower  bound  on  the  variance  of  an  unbiased 
statistical  estimator.  The  reciprocal  of  the  Fisher  Information,  J{x',6),  of  the  parameter  0,  is  the 
lower  bound  on  the  variance  of  an  unbiased  estimator  of  the  parameter  0. 

cov{^}>j‘  (B.l) 

Consider  the  problem  of  estimating  a  parameter  vector  0  from  sample  values  of  an  n-variate 
time  series  x{t) .  The  elements  of  the  Fisher  Information  for  an  unbiased  estimator  0  are  given 
by: 

J,{x-e)  =  -E 

Where  p{x;d)  is  the  conditional  probability  density  function.  For  the  case  where  the  elements  of 
X  are  jointly  Gaussian,  Bangs  has  shown  that  the  elements  of  the  Fisher  Information  matrix  are 
given  by: 


Where  R  is  the  data  covariance  matrix  given  hy  R  =  E {xx^}  .  The  FI  evaluation  in  (B.3)  can  be 
simplified  if  the  matrix  R  is  diagonal.  This  is  achieved  by  replacing  the  data  vector  x  by  its 
Fourier  coefficients,  which  are  asymptotically  uncorrelated. 

Let  x{t)  be  a  discrete  time  stationary  zero-mean  n-variate  process  with  autocorrelation  function 
R(t)  =  E {x{t)x^{t  +  r)} .  Its  power  spectrum  density  is  defined  as: 


\^R-^^rA  (B.3) 

^0,  de, 


d^\Ap{x\e)  ^ 
dOBO.,  I  ^  ^ 

*  J  J 


Sico)=  (B.4) 

r=— 00 


The  elements  of  the  Fisher  Information  from  (B.3)  can  then  be  expressed  as: 


Jij  = - J  tr 


do, 

V  I 


(B.5) 
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This  result  also  holds  for  a  band-limited  continuous  time  process  which  is  sampled  at  the  Nyquist 
rate. 

For  the  case  of  time  delay  estimation,  consider  two  noisy  signals  r^{t)  and  r^{t)  which  contain  a 
common  stationary  process,  s{t),  delayed  in  time. 

r,{f)  =  s{t)  +  nAt) 

'  ‘  (B.6) 

r2(t)  =  s{t-D)  +  ^2(0 

where  n^{t)  and  ^2(0  ^re  uncorrelated  additive  noise  processes.  The  relative  delay,  D,  is 
estimated  from  measurements  {r^{t),r2{t),\  <t<T).  The  Cramer-Rao  bound  for  delay  estimation 

error  is  formulated  by  first  forming  the  spectral  matrix  of  the  vector  process  . 


S{(o)  +  N,{(o)  S(co)e^‘“^ 


(B.7) 


Then,  upon  differentiating  the  above  expression  with  respect  to  the  time  delay  parameter  and 
performing  the  necessary  matrix  multiplications,  the  following  is  found; 


a)  a)  |S(®)|L  0  q/S\o)\ 


(B.8) 


Where  |S((y)|  is  the  determinant  of  the  spectral  matrix  S(o).  When  (B.8)  is  inserted  into 
equation  (B.5),  the  resulting  expression  for  the  Fisher  Information  is: 


T 


2k 


^^S(a)l 


N2((a)J 


I  -I- 


,)  + 


-dco 


N2(a>)J 


(B.9) 


From  equation  (B.9)  we  can  see  that  the  quadratic  weighting  in  the  Fisher  Information  arises  as  a 
result  of  differentiating  S{a))  with  respect  to  the  time  delay  parameter,  D. 


B.2  The  effects  of  carrier  phase  synchronization  on  the  estimation  process 


This  appendix  demonstrates  the  effects  of  synchronization  errors  (between  the  local  oscillator 
and  the  received  carrier)  on  the  TDOA/FDOA  estimation  process.  Given  two  complex  baseband 
signals  with  a  common  component  shifted  in  frequency  and  time: 


Xj(t)  = 


(B.IO) 


The  TDOA/FDOA  estimation  process  is  defined  as: 
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max  I  x^{t)x*2{t+  T)e  (B.ll) 

or  equivalently: 

max  }  s{t)e^^’^^‘s{t+  (B.12) 

0 


Assuming  that  the  signal  is  demodulated  without  the  use  of  earrier  phase  synchronization,  then 
let  the  phase  offset  at  sensor  i  be  6^ .  The  baseband  representations  of  the  received  signals  at  the 
two  sensors  are  now: 


(B.13) 

The  TDOA/FDOA  estimation  process  requires  searching  over  the  magnitude  of  the  modulated 
cross-correlation: 

1 

0 

(B.14) 

0 

(B.15) 

0 

(B.16) 

=  \  (B.17) 

0 

Therefore,  TDOA/FDOA  parameter  estimation  does  not  require  the  use  of  a  coherent  phase 
reference  when  removing  the  carrier  component  from  the  observed  signal.  While  an  accurate 

phase  reference  is  necessary  for  symbol  detection,  it  has  no  effect  on  the  TDOA/FDOA 
estimation  process. 


53 


Appendix  C:  Experimental  Setup 

C.l  Preliminaries 


The  experiments  performed  in  this  thesis  require  the  use  of  the  MATLAB  eore,  Communications 
Toolbox,  and  Signal  Processing  Toolbox.  The  following  m-files  were  used  in  the  simulations. 

crMsg.m,  tdoa.m,  codec. m,  bitalloc.m,  bitallocfi.m,  addnoise.m,  Tx  symbols.m 

C.l  Decimation  Simulations 


This  simulation  decimates  and  interpolates  digital  linearly  modulated  signals  of  different  symbol 
rates,  generating  100  random  messages  for  each  specified  symbol  rate.  Decimation  and 
subsequent  interpolation  are  performed  for  multiple  factors.  The  TDOA  parameter  is  then 
estimated  using  both  the  compressed  and  uncompressed  messages.  The  results  for  each  message 
will  be  saved  in  the  directory  specified  by  dec _path. 

%test  effects  of  downsampling  (varying  symbol  rate  and  decimation  factor) 
%Preliminaries : 

%  msgR  -  message  parameter  structure 
%  P,u  -  sensor/emitter  coordinates 

%  cdcP  -  codec  parameter  structure 

%Experimental  results: 

%  td  -  matrix  containing  TDOA  estimates  obtained  from  uncompressed  signals 
%  ctd  -  TDOA  estimates  obtained  from  compressed  signals 


dec_path= ' results/dec/ '  ; 

M=[2:2:32  35  36  40  42  45  48]; 

F=[2  2  3:5  4  7  4  6  5  11  6  13  7  6  8 
F=[F;  12222324342424 
N  sym=256 : 256 : 12800; 
R“sym=(50:50:2500)  *1000; 

for  e=l: length (R  sym)  , 
msgR.R  sym=R  sym(e); 
msgR.N  sym=N  sym(e); 
msgR. x= [ ] ; 


%path  in  which  results  are  saved 
%interpolation  factors 
768798];  %decimation  factors 

54565656];  %(two  stages) 
%number  of  symbols 
%symbol  rates 

%for  each  symbol  rate 

%clear  symbol  vector 


clear  rtd  ctd  td 
for  k=l:100, 

[msgk  msgP] =crMsg (u, P, msgR, [ ] ) ; 


%repeat  100  times 
%create  new  message 


msgkc=msgk; 


for 


j=l: 

for 


: length (M) , 

i=l: (size(msgk,l)-l), 

K=length (msgk{ i, 1 } ) ; 

Pad=ceil (K/M(j))*M(j)-K; 
x= [msgk{ i, 1 } ;  zeros ( Pad, 1 ) ] 

%decimate  and  interpolate 
msgkc { i , 1 } =interp (decimate (decimate (x,F(l,j) , 'FIR 
msgkc { i , 1 } =msgkc { i , 1 } ( 1 : end-Pad) ;  %remove 


%copy  other  stuff  over 
%for  each  decimation  factor 


jcalculate 

%zero-pad 


pad 

signal 


)  ,F(2,  j) 
pad 


'FIR' ) ,M(j) ) 
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end 

ctd  ( j  ,  : ,  k)  =tdoa  (msgk;c,msgP,  1)  ;  %compressed  TDOA  estimate 
[td(j,:,k)  rtd ( j  ,  ; , k) ] =tdoa (msgk, msgP, 1 ) ;  %uncompressed  TDOA  estimate 
%display  status 

disp ( [num2str (M ( j ) )  '  num2str (ctd ( j , : , k) -td ( j , : , k) ) ] ) ; 

end 

save ( [dec_path  num2str (e) ] , ' ctd ' , ' td ' , ' M ' , ' msgP ' ) ;  %save  results 

end 

end 


C.3  Subband  Encoding  Simulations 


This  simulation  will  generate  100  random  messages  and  perform  subband  eneoding  using  MSB 
and  FI-based  distortion  criteria.  A  bit  allocation  will  be  computed  for  each  intercepted  signal. 
The  TDOA  parameter  is  estimated  using  signals  produced  by  both  compression  methods.  The 
results  are  saved  to  the  fde  sub  Jile. 


%test  FI  vs  MSE  bitalloc  with  additive  noise  (varying  SNR) 

%Preliminaries : 

%  msgR  -  message  parameter  structure 
%  P,u  -  sensor/emitter  coordinates 

%  cdcP  -  codec  parameter  structure 

%Experimental  results: 

%  td  -  matrix  containing  TDOA  estimates  obtained  from  uncompressed  signals 
%  ftd  -  TDOA  estimates  for  signals  compressed  using  FI  criteria 

%  mtd  -  TDOA  estimates  for  signals  compressed  using  MSE  criteria 

%  bf  -  subband  bit  allocations  using  FI  criteria 

%  bm  -  subband  bit  allocations  using  MSE  criteria 


clear  td  ftd  mtd  bm  bf 

sub  f ile= ' results/ subband/ experiment! ' ; 

N=fliplr( [0:5:40] ) ;  %SNR 

for  i=l:100,  %repeat  100  times 

[msgb  msgP] =crMsg (u, P, msgR, [ ] ) ;  %create  new  message 

for  j =1 : length (N) ,  %for  each  SNR 

cdcP.n=max  (abs  (msgb{  1,  1 } )  )  "'2/10'^  (N  ( j  ) /lO)  ;  %noise  variance 
msgk=addnoise (msgb, N ( j ) , N ( j ) ) ;  %add  noise  (WON) 

cdcP . op= ' MSE ' ;  %change  mode  to  MSE 

[msgkc  bm ( i , j ,:,:)] =codec (msgk, cdcP) ;  %compress  using  MSE  alloc 
msgkc { 3 , 1 } =msgk{ 3 , 1 } ;  %dont  compress  signal  3 

mtd (i, :, j ) =tdoa (msgb, msgP, 1) -tdoa (msgkc, msgP, 1) ;  %MSE  TDOA  estimates 
cdcP . op= ' FI ' ;  %change  mode  to  FI 

[msgkc  bf ( i , j ,:,:)] =codec (msgk, cdcP) ;  %compress  using  FI  alloc 
msgkc { 3 , 1 } =msgk{ 3 , 1 } ;  %dont  compress  signal  3 

ftd (i, :, j ) =tdoa (msgb, msgP, 1) -tdoa (msgkc, msgP, 1) ;  %FI  TDOA  estimates 


td ( i ,:, j ) =tdoa (msgb, msgP, 1 ) -tdoa (msgk, msgP, 1 ) ;  %no-comp  TDOA  estimates 
disp(['MSE:  '  num2str (mtd (i,  :  ,  j  )  )  '  FI:  '  num2str (ftd (i,  :  ,  j  )  )  '  NO:  ' 

num2str (td (i, : , j ) ) ] ) ;  %display  status 

end 

save ( sub_f ile,  ' f td '  ,  ' mtd ' ,  ' td '  ,  ' msgP  '  ,  ' bm ' ,  ' bf ' ) ;  %save  results 

end 
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This  simulation  will  generate  100  random  messages  and  perform  subband  encoding  using  MSB 
and  FI-based  distortion  criteria;  however,  unlike  in  the  previous  experiment,  the  bit  allocation 
will  be  computed  only  once  using  a  randomly  generated  message.  This  simulation  generates  the 
results  for  the  suboptimal  approach  that  was  presented  in  section  4.3. 

%test  FI  vs  MSE  bitalloc  with  additive  noise  (varying  SNR,  one-time  alloc) 

%  msgR  -  message  parameter  structure 
%  P,u  -  sensor/emitter  coordinates 

%  cdcP  -  codec  parameter  structure 

%Experimental  results: 

%  td  -  matrix  containing  TDOA  estimates  obtained  from  uncompressed  signals 
%  ftd  -  TDOA  estimates  for  signals  compressed  using  FI  criteria 

%  mtd  -  TDOA  estimates  for  signals  compressed  using  MSE  criteria 

%  bf  -  subband  bit  allocations  using  FI  criteria 

%  bm  -  subband  bit  allocations  using  MSE  criteria 

clear  td  ftd  mtd 

sub  f ile= ' results/ subband/ experiment! ' ; 

N=fliplr( [0:5:40] ) ;  %SNR 


[msgb  msgP] =crMsg (u, P, msgR, []) 


for 


j  =1 : length (N) , 

cdcP.n=max (abs (msgb{ 1, 1 } ) ) ^2/10 
msgk=addnoise (msgb, N ( j ) , N ( j ) ) ; 
cdcP . op= ' MSE ' ; 

[msgkc  bm ( j , : , : ) ] =codec (msgk, cdcP) 
cdcP . op= ' FI ' ; 

[msgkc  bf ( j , : , : ) ] =codec (msgk, cdcP) 


%create  message 
%for  each  SNR 
(N(j)/10);  %noise  variance 
%add  noise  (WON) 


for  bit  alloc 


%change  mode  to 
%compress  using 
%change  mode  to 
%compress  using 


MSE 

MSE  alloc 
FI 

FI  alloc 


end 


cdcP . op= [ ] ; 
for  i=l:100, 

[msgb  msgP] =crMsg (u, P, msgR, []) 
for  j =1 : length (N) , 

msgk=addnoise (msgb, N ( j ) , N  ( 
cdcP . b=squeeze (bm ( j , 1 , : ) ) ; 
msgkc=codec (msgk, cdcP) ; 
msgkc { 3 , 1 } =msgk {3,1}; 
mtd ( i , : , j ) =tdoa (msgb, msgP, 
cdcP . b=squeeze (bf ( j , 1 ,  : ) )  ; 
msgkc=codec (msgk, cdcP) ; 
msgkc { 3 , 1 } =msgk { 3,1}; 
ftd(i, :, j ) =tdoa (msgb, msgP, 
td ( i , : , j ) =tdoa (msgb, msgP, 1 
disp(['MSE:  '  num2str(mtd( 
num2str (td (i,  :  ,  j )  )  ]  )  ; 
end 

save (sub_file, ' ftd' , 'mtd' , ' td' 
end 


%repeat  100  times 
;  %create  new  message 

%for  each  SNR 

j));  %add  noise  (WON) 

%load  MSE  bit  allocation 
%compress  using  MSE  alloc 
%dont  compress  signal  3 
1 ) -tdoa (msgkc, msgP, 1 ) ;  %MSE  TDOA  estimates 
%load  FI  bit  allocation 
%compress  using  FI  alloc 
%dont  compress  signal  3 
1 ) -tdoa (msgkc, msgP, 1 ) ;  %FI  TDOA  estimates 
) -tdoa (msgk, msgP, 1 ) ;  %no-comp  TDOA  estimates 
i,:,j))  '  FI:  '  num2str (ftd (i,  :  ,  j  )  )  '  NO:  ' 

%display  status 


, ' msgP ' , ' bm ' , ' bf ' ) ; 


%save  results 
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